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Abstract 

The understanding of the structural and thermal properties of membranes, low-dimensional flexi- 
ble systems in a space of higher dimension, is pursued in many fields from string theory to chemistry 
and biology. The case of a two dimensionai (2D) membrane in three dimensions is the relevant one 
for dealing with real materials. Traditionally, membranes are primarily discussed in the context of 
biological membranes and soft matter in general. The complexity of these systems hindered a real- 
istic description of their interatomic structures based on a truly microscopie approach. Therefore 
theories of membranes were developed mostly within phenomenological models. From the point of 
view of statistical mechanics, membranes at finite temperature are systems governed by interacting 
long-range fluctuations. 

Graphene, the first truly two-dimensional system consisting of just one layer of carbon atoms, 
provides a model system for the development of a microscopie description of membranes. In the 
same way that geneticists have used Drosophila as a gateway to probe more complex questions, 
theoretical chemists and physicists can use graphene as a simple model membrane to study both 
phenomenological theories and experiments. In this Account, we review key results in the mi- 
croscopie theory of structural and thermal properties of graphene and compare them with the 
predictions of phenomenological theories. The two approaches are in good agreement for the var- 
ious scaling properties of correlation functions of atomic displacements. However, some other 
properties, such as the temperature dependence of the bending rigidity, cannot be understood 
based on phenomenological approaches. We also consider graphene at very high temperature and 
compare the results with existing models for two-dimensional melting. The melting of graphene 
presents a different scenario, and we describe that process as the decomposition of the graphene 
layer into entangled carbon chains. 



I. INTRODUCTION 



Understanding the structural and thermal properties of two dimensionai (2D) systems 
is of great interest in many fields including mechanics, statistical physics, chemistry and 
biologyR Traditionally, it was discussed mainly in the context of biological membranes and 
soft condensed matter. The complexity of these systems hindered any truly microscopie 
approach based on a realistic description of interatomic interactions. Phenomeno logicai 
theories of membranes^ based on elasticità reveal nontrivial scaling behavior of physical 
properties, like in- and out-of-plane atomic displacements. In three-dimensional (3D) sys- 
tems, this type of behavior takes place only close to criticai points^, whereas in 2D this occurs 
at any finite temperature. The discovery of graphene^, the first truly 2D crystal made of 
just one layer of carbon atoms, provides a model system for which an atomistic description 
becomes possible. The interest for graphene has been triggered by its exceptional electronic 
properties (for review see^) but the experimental observation of ripples in freely suspended 
graphene^ has initiated a theoretical interest also in the structural properties^. Ripples 
or bending fluctuations have been proposed as one of the dominant scattering mechanisms 
that determine the electron mobility in graphene^. Last but not least, the structural state 
influences the mechanical properties that are important for numerous potential applications 
of graphene^^. 

Graphene is a crystalline membrane, with finite resistance to in piane shear deformations, 
contrary to liquid membranes as soap films. Moreover, for graphene, this resistance is 
extremely high since the carbon-carbon bond is one of the strongest chemical bonds in 
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nature. The Young modulus per layer of graphene is 350 N/m, ari order of magnitude larger 
than that of stee l 10 * 11 *. Phenomenological theories just assume that the membrane thickness 
is negligible in comparison to the lateral dimensions. Once a microscopie treatment is 
allowed, one can also distinguish between, e.g., single layers and bilayers^. 

The aim of this review is to summarize the contribution of microscopie treatments of 
graphene as the simplest (prototype) membrane to our general understanding of 2D systems. 
To this purpose we first review the main results of phenomenological theories, particularly 
those that can be directly compared to results of atomistic approaches. 

II. PHENOMENOLOGICAL THEORY OF CRYSTALLINE MEMBRANES 

The standard theory of lattice dynamics is based on the harmonic approximation as- 
suming atomic displacements from equilibrium to be much smaller than the interatomic 
distance d. For 3D crystals, this assumption holds up to melting according to the empirical 
Lindemann criterion. For 2D crystals the situation is so different that Landau and Peierls 
suggested in the 1930's that 2D crystals cannot exist. Later, their qualitative arguments 
were made more rigorous in the context of the so-called Mermin- Wagner theorem (see ref- 
erences irP. Since graphene is generally considered to be a 2D crystal this point needs to be 
clarified first. 

A. Lattice dynamics of graphene 

By aauming that atomic displacements u satisfy the condition 

< <i >« d 2 (1) 

where n labels the elementary celi and j the atoms within the elementary celi, we can expand 
the potential energy V(R) up to quadratic terms (harmonic approximation): 

v(RrJ)=v{^)+\ £ < n ,,<f<? (2) 

n,n' i,j ufi 

where the matrix A is the force Constant matrix, R n j = + u n j and R^j = f n + pj where 
f n are the vectors of the 2D Bravais lattice and pj are basis vectors. Lattice vibrations 
are then described as superposition of independent modes, called phonons, characterized by 
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wavevector q and brandi number £ = 1, 3z/ where v is the number of atoms per unit cell. 
The squared phonon frequencies oo 2 {q) are the eigenvalues of the 3u x 3u dynamical matrix 



A afì 

(3) 

n V 1 J 

where Mj is the mass of atom j. For graphene, Mj = M is the mass of the carbon atom 
and by symmetry AfJ = = and = D^ 2 - Translational invariance requires that no 
forces result from a rigid shift of the crystal, implying: 

£<* = ° ( 4 ) 

whence 

D#(?=0) 0) = (5) 
Therefore, there are six phonon branches in graphene: 

1) The acoustic flexural mode ZA 

" 2 ZA {q) = D{{{q) + DH{q) (6) 

2) The optical flexural mode ZO (tt||Oz) 

u 2 zo{q) = Dll{q)-DH{q). (7) 

3), 4) Two acoustic in-plane modes, with u 2 (q) equal to the eigenvalues of the 2x2 matrix 

Dtf(fì+ng® (a,P = x,y) (8) 

5), 6) Two optical in-plane modes, with co 2 (q) equal to eigenvalues of the 2x2 matrix 

Dt (q) - Z#(g) (a,P = x,y) (9) 

If the 2D wavevector q lies in symmetric directions, branches (3)- (6) can be divided into 
longitudinal e\\q and transverse e _L q modes. Due to [H] for acoustic modes oo 2 oc q 2 at 
q — > 0. For the ZA mode, however, the terms in q 2 disappear as well and oozaÌq) 00 
This follows from the invariance with respect to rotations of a 2D crystal as a whole in the 
3D space, namely for uniform rotations of the type 

u nj = ò<\>m x R^, (10) 



FIG. 1. Phonon spectrum of graphene calculated with LCBOPII. Adapted fromP^l 




where 5(p is the rotation angle and rh the rotation axis in the xy-p\ane. These rotations 
should not lead to any forces or torques acting on the atoms. Hence, 

J2 A °Ufnrn = = *,!/)■ (H) 

It follows from 11 3] that 



g^lDSW + BSOA^O d2) 
and, thus, the expansion of[6]starts with terms of the order of g 4 ; therefore, 

Uza oc q 2 (13) 

at q — > 0. The very low frequency of uza for q — > has important consequences for the 
stability and thermal properties as we discuss next. In[T]we show the phonon spectrumP^l 
calculated with the so-called long-range carbon bond order potential (LCBOPII)P^ used in 
the atomistic simulations presented later. 

Let us consider now the case of finite temperatures. In the harmonic approximation, the 
mean-square atomic displacement is 

<^>=Eto(«&)*K)-*(^) <"> 

where A = (q,£) are phonon labels, e is the polarization vector and No is the number of 
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elementary cells. For in-plane deformations at any finite temperature the sum in 14 



logarithmically divergent due to the contribution of acoustic branches with oj oc q for q — > 0. 
This divergence is cut at the minimal wavevector q min ~ L^ 1 (L is the sample size), thus 

< > = < Vi, >« ^5 !" (§) (15) 

where c s is the average sound velocity. This result led Landau and Peierls to the conclusion 
that 2D crystals cannot exist. Strictly speaking, this means just the inapplicability of the 
harmonic approximation, due to violation of [T] A more rigorous treatment, however, does 
confirm this conclusion (see 2 ). For a = z, the situation is even worse, due to the much 
stronger divergence of ZA phonons ([6]). One can see from 14 that 



< hi, >oc 



where E at is of the order of the cohesive energy. Henceforth we use the notation h = u z 



and denote u = (u x ,u y ) as a 2D vector. 



B. The statistical mechanics of crystalline membranes 

We have shown that the harmonic approximation cannot be applied at any finite tem- 
perature to 2D crystals neither for in-plane nor for out-of-plane modes since the condition 
[T] is violated due to divergent contributions of acoustic long-wavelengths modes with q — > 0. 
In this situation, it becomes necessary to consider anharmonic interactions between in-plane 
and out-of-plane modes. In the limit q — > 0, acoustic modes can be described by elasticity^. 
The corresponding effective Hamiltonian % reads 

n= l -jd 2 x( K K {v 2 hf + + , (ì?) 

where the deformation tensor u a p is 

1 / dup du a dh dh\ 
al3 2 \dx a dxp dxadxp) 

k is the bending rigidity and fi and A are Lamé coefficients. In the deformation tensor, we 

have kept the nonlinear terms in dh/dx a but not du^/dx a since out-of-plane fluctuations 



are stronger than in-plane ones (compare 15 and 16). If we neglect ali nonlinear terms in 



the deformation tensor, then % in q representation becomes: 

Ho = ^E g%f + \ E W\*# + ( A + /*) • ^) 2 ] ( 19 ) 

q g 



where the subscript indicates the harmonic approximation and h$ and u$ are Fourier 
components of h(f) and u(r), respectively, with f= (x,y). 
The correlation functions in harmonic approximation are 

G (q) =< \h g f > = (20) 



Df(q) =< u* a *u M > 



T 



q 2 (A + 2/i)g 2 



QaQp 



(21) 



where <>o means averaging with the Hamiltonian T-Lq (19). 



For a surface z = h(x,y) the components of the normal are 

dh 1 



n. 



n. 



n. 



dx y/i + |V/i| 2 
dh 1 



dy^/T+WW 
1 



(22) 
(23) 
(24) 



where Vh is a 2D gradient. If \Vh\ << 1, the normal-normal correlation function is related 



tO < |/ig-| 2 > 



(25) 



On substituting 20 into 25 we find 



< ntfi-ff > = 



T 

Kq 2 



(26) 



A membrane is globally fiat if the correlation function < n^n^ > tends to a Constant as 



R — > oo (normals at large distances have, on average, the same direction). 26, instead, 
leads to a logarithmic divergence of < n Q n^ >. Moreover, the mean square in-plane and 



out-of-plane displacements calculated from 20 and 21 are divergent as L — > oo as already 
shown. Again, we conclude that the statistical mechanics of 2D systems cannot be based 
on the harmonic approximation. Taking into account the coupling between u and h due to 



the nonlinear terms in the deformation tensor 18 drastically changes this situation. We can 
introduce the renormalized bending rigidity KR,(q) by writing 

T 



G{q) 



KR(q)q 4 



(27) 



The fìrst-order anharmonic correction to k is 



3TY 

5k = K R (q) - k = = (28) 



where Y = \ { ^ ] is the 2D Young modulmPl At 



3TY , . 

V 07TK Z 

the correction Sk = k, and the coupling between in-plane and out-of-plane distortions cannot 
be considered as a perturbation. The value q* plays the sanie role as the Ginzburg criteriorP 
in the theory of criticai phenomena: below q* interactions between fluctuations dominate. 
Note that in the theory of liquid membranes, there is also a divergent anharmonic correction 
to k of completely different origirP 

* (il (3°) 



Ah \qd 



This term has sign opposite to the one of a crystalline membrane, [28} and is much smaller 
than the latter. 

In presence of strongly interacting long-wavelength fluctuations, scaling considerations 
are extremely usefuP. Let us assume that the behavior of the renormalized bending rigidity 
Kr{q) a t small q is determined by some exponent 77, Kn{q) oc q~ v yielding 



A A 
q 4 »»g ' q 2 % 



where the parameter q = \/Y / k of the order of d~ l is introduced to make A dimensionless. 
One can assume also a renormalization of the effective Lamé coefficients A^(g), ^n{q) oc q Vu 
which means 

<<^>oc-^- (32) 



Finally, we assume that anharmonicities change 16 into 



< h 2 >oc L 2C - (33) 

The values r],T] u and ( are similar to criticai exponents in the theory of criticai phenomena. 
They are not independenlP^ 



C = 1 - 77/2, Vu = 2 - 277 
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(34) 



The exponent r\ u is positive if < r] < 1. The so-called Self-Consistent-Screening- 
Approximatiori^ gives rj ~ 0.82 whereas a more accurate renormalization group approacfP^ 
yields r\ rj 0.85. This means that, interactions make out-of-plane phonons harder and 
in-plane phonons softer. 



that 



The temperature dependence of the Constant A in [27] can be found from the assumption 
should match at q = q*, giving A = a (T/k) where a is a dimensionless 
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and 



31 



factor of the order of one. 

Now we are ready to discuss the possibility of long-range crystal order in 2D systems at 
finite temperatures. The true manifestation of long-range order is the existence of delta- 
function (Bragg) peaks in diffraction experiments. The scattering intensity is proportional 
to the structure factor 



S ÌQ) = ^2 ( eXP wi^ni ~ Èn'j'j ) ( 35 ) 



nn' jj' 

that can be rewritten as 

S (0) = ^ eX P W (*n - T n ')] 22 eX P W(Pj ~ P?)\ W nj,n'j' (36) 
nn' jj' 

where 

W njtn 'f = (exp [iq(u nj - u n >f))} (37) 

In 3D crystals, one can assume that the displacements u n j and u n 'f are not correlated for 
\r n — f n '\ — > oo so that 

W nj,n'j' = (exp (iqunj)) (exp (-iqu n/j/ )) = rrij{q)m*,{q) (38) 

where rrij(q) are Debye-Waller factors that are independent of n due to translational invari- 
ance. Therefore, for q = g (reciprocai lattice vectors), where exp (iqr n ) = 1, the contribution 
to S(q) is proportional to Nq , whereas for a generic q it is of the order of A^ . The Bragg 
peaks at q = g are, therefore, sharp; thermal fluctuations decrease their intensity (by the 
Debye-Waller factor) but do not broaden the peaks. The observation of such peaks is an 
experimental manifestation of long-range crystal order. In 2D the correlation functions of 
atomic displacements do not vanish as \r n — f n >\ — > oo. Indeed, in the continuum limit, 
u n j — » (u(r), h{f)) and we have 

(\h(f) - /i(r)f) =2^(|/i(gì 2 ) [1 - cos (<f (f - r'))] ~ |f - 7\ 2<: (39) 



1 2 



u(r) 



u(r 



2 (\u(q\ 2 ^ [1 — cos (<f (f — f ))] ~ |r — r , |' 7 " 



(40) 



after substitutions of 3~ì~j32f . Thus the approximation 38 does not apply. As a result, the 
sum over n' at a given n is convergent, and S(q = g) oc A^; instead of a delta-function Bragg 
peak we have a sharp maximum of finite width. This means that, rigorously speaking, the 
statement that 2D crystals cannot exist at finite temperatures is correct. However, the 
structure factor of graphene stili has sharp maxima at q = g and the crystal lattice can 
be determined from the positions of these maxima. In this restricted sense, 2D crystals do 
exist, and graphene is a prototype example of them. 

It was found experimentally by transmission electron microscopy, that freely suspended 
graphene at room temperature is rippled 6 . The existence of these thermally induced ripples 
motivated our atomistic Monte Carlo simulations^^ summarized in the next section. 



III. ATOMISTIC SIMULATIONS OF STRUCTURAL AND THERMAL PROP- 
ERTIES OF GRAPHENE 

As discussed before the thermal properties of 2D crystals are determined by longwave- 
length fluctuations. Therefore, one needs to deal with large enough systems to probe the 
interesting regime of strongly interacting fluctuations. This requirement rules out, in prac- 
tice, first principle approaches in favor of accurate empirical potentials. The unusual struc- 
tural aspects of graphene, make it desirable to describe different structural and bonding 
configurations, beyond the harmonic approximation, by means of a unique interatomic po- 
tential. Bond order potentials are a class of empirical interatomic potentials designed for 
this purpose (see 15 and references therein). They aim at describing also anharmonic effects 
and the possible breaking and formation of bonds in structural phase transitions. They 
allow to study without further adjustment of parameters, ali carbon structures, including 
the effect of defects, edges and other structural changes, also as a function of temperature 
as well as phonon spectra. We have used the so-called long-range carbon bond order poten- 
tial LCBOPII^. Its main innovative feature is the treatment of interplanar van der Waals 
interactions, that allows to deal with graphitic structures. To calculate equilibrium prop- 
erties as a function of temperature, we have performed Monte Carlo simulations either at 
Constant volume or Constant pressure. We first discuss the results for correlations functions 

10 



that can be directly compared to the scaling behavior discussed previously. Then, we report 
the temperature dependence of several structural properties of graphene. Lastly we discuss 
the melting of graphene in relation to its 3D counterpart, graphite, and to 2D models of 
melting. 



A. Structural properties and scaling 

We compare the results of atomistic Monte Carlo simulations to the scaling behavior of 



G(q) (31 ). From 31 , one can see that G(q) can be calculated in two ways, either by calculating 
directly the correlator of h(f) or of the normal n(f). In doing this, it is important to have 
h(r) and n(f) calculated at lattice sites smoothened by averaging over the neighbors, as 



described in detail in RefP^. Only by such a procedure one verifies numerically 25 forg < 10 



nm 1 which gives the limit of applicability of a continuum description to graphene. The 



interesting regime is q << q* (29). For graphene at room temperature q* = 2.4 nm -1 . Since 
simulations are done for samples of dimension L x x L y with periodic boundary conditions, 
the smallest values of q that can be reached are 27r/L x and 2ir /L y . For the largest samples, 
we have found that straightforward Monte Carlo simulations based on individuai atomic 
moves, could not provide enough sampling for the smallest wavevectors. For this reason in 
our first paper on ripples in graphene^ we were not able to check the scaling laws in the 
anharmonic regime. Later, we have reached this regime by devising a numerical technique 
that we have called wave moves, where collective sinusoidal long wavelengths displacements 
of ali atoms where added in the Monte Carlo equilibration procedure^. In [2] we present 
< \rìq\ 2 >= q 2 G(q) calculated with wave moves which displays a clear change of the slope 
\n(G(q)) versus ln(g) around q ~ q*. Notice that < \n^\ 2 > grows for q > 10 nm -1 reaching 
a maximum at the first Bragg peak. According to the phenomenological theory described 
before, the change of scaling behaviour at q << q* is related to the coupling of in-plane and 



out-of piane fluctuations. To check this, we also show in [3] the correlation function T(q) 

r{q)=<(u x ) q ih 2 U> (41) 

which becomes almost zero at q > q*. For smaller samples the coupling is reduced as 
expected for a property that is determined by the region of long wavelengths fluctuations. 
The temperature dependence of the bending rigidity k(T) can be extracted from the 
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FIG. 2. Normal-normal correlation function q 2 G(q) for three samples with indicateci number of 
atoms N. For the largest, N = 37888, L x = 314.82 À, L y = 315.24 À. Adapted from RefP 
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FIG. 3. The function T(q) for three values of N 
< |n<?| 2 > using 26 The resulti are shown in|4|where one can see the rapid growth with 



temperature. This effect should not be confused with the correction 5k of 28 since the latter 
is strongly g-dependent. The temperature dependence of k of|4]cannot be described within 
the Self-Consistent Screening Approximation for our model Hamiltonian [ì~7 { 9 . 

The temperature dependence of «, as that of ali parameters of phonon spectraP^, is an 



anharmonic effect that goes beyond the model 17, namely it does not result from the coupling 



of acoustic out-of piane phonons with acoustic in-plane phonons only. Other anharmonicities, 
like coupling to other phonons have to be invoked. This is an example of effects that can 
be studied within atomistic simulations but not within the elasticity theory. The potential 
energy given by LCBOPII includes by construction anharmonic effects. 
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FIG. 4. Temperature dependence of k as found by fitting < {nA 2 > to 31 Adapted from Re# 



FIG. 5. Temperature dependence of the lattice parameter a calculated by Monte Carlo simulations 
at zero pressure. Adapted from RefP 
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Purely anharmonic effects are the temperature dependence of lattice parameter and elastic 
moduli. The temperature dependence^ of the lattice parameter a is shown in [H] and those 
of the shear modulus /1 and adiabatic bulk modulus b A = X + fi in [6j The most noticeable 
feature in [5] is the change of sign of da/dT, namely a change from thermal contraction to 
thermal expansion around 1000 K. 

Usually thermal expansion is described in the quasiharmonic approximatiorp3 where the 
free energy is written as in the harmonic approximation but with volume dependent phonon 
frequencies u\. This dependence is described by the Grùneisen parameters 

Sin u\ 



7a 



(42) 
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FIG. 6. Temperature dependence of the bulk modulus Òa and shear modulus //. Adapted from 
RefP 

where Q is the volume (the area for 2D systems). In most solids, phonon frequencies grow 
under compression, which corresponds to positive Grùneisen parameters and thermal ex- 
pansion. Graphene and graphite are however exceptional, as illustrated in [7] presenting the 
corresponding calculation with LCBOPlP^. One can see that both the ZA and ZO branches 
have 7 < almost in the whole Brillouin zone as found already, within density functional 
calculations^. 

Experimentally, graphite has a negative thermal expansion coefficient up to 700KP. This 
behavior has been explained in the quasiharmonic approximation in Ref.^. For graphene, 
they predicted negative da/dT at ali temperatures. Negative thermal expansion of graphene 
at room temperature has been confirmed experimentally^. The linear thermal expansion 
coefficient was about -1CT 5 K _1 , a very large negative value. According to the quasiharmonic 
theory, it was found to be more or less Constant up to temperatures of the order of at least 
2000 K ,in contrast to the atomistic simulations of [5] Thus the change of sign of da/dT 
should be attributed to self-anharmonic effects^, namely to direct effects of phonon-phonon 
interactions. Very recently, it was confirmed experimentally that da/dT, while remaining 
negative, decreases in modulus with increasing temperature up to 400K 131 , which can be 
considered as a partial confirmation of our prediction. 

Also the temperature dependence of the shear modulus fi, shown in [6] is anomalous 
since typically dfi/dT < at any temperature. The change of sign of dfi/dT < occurs 
roughly at the same temperature of da/dT. The room-temperature values of the elastic 

14 



2 




Wavevector 




FIG. 7. Top) Grùneisen parameters calculated for graphene with LCBOPI1 14 ; bottoni) Phonon 
spectrum of graphene for the equilibrium value of the interatomic distance 1.42 A (red solid), and 
two larger values, 1.43 A (blue dashed) and 1.4^ À (green dotted). Courtesy of L.J. Karssemeijer 



constants are /i ~ 10 eV/À 2 and Òa ~ 12 eV/A 2 . The corresponding Young modulus Y lies 
within the error bars of the experimental value^ Y ~ 340 ± 50 Nm _1 . The Poisson ratio 
v — ipA ~ aO/(&a + A*) is found to be very small, of the order of 0.1. 

B. Melting of graphene 

Melting in 2D is usually described in terms of creation of topological defects, like un- 
bound disclinations that destroy orientational order and unbound dislocations that destroy 
translational order^. In the hexagonal lattice of graphene, typical disclinations are pen- 
tagons (5) and heptagons (7) while dislocations are 5-7 pairs. Our atomistic simulations^ 
have given an unexpected scenario of the melting of graphene as the decomposition of the 
2D crystal in a 3D network of 1D chains. A cruciai role in the melting process is played 
by the Stone- Wales (SW) defects, non-topological defects with a 5-7-7-5 configuration. The 
SW defects have the smallest formation energy and start appearing spontaneously at about 
4200 K. It is the clustering of SW defects that triggers the spontaneous melting around 4900 
K in our simulations. 

In[8]we show a typical configuration of graphene on the way to melting at 5000 K. The 
coexistence of crystalline and molten regions indicates a first order phase transition. The 
most noticeable features are the puddles of graphene that have molten into chains. The 
molten areas are surrounded by disordered 5-7 clusters, resulting from the clustering and 
distortion of SW defects. Isolated and pairs of SW defects are also present whereas we 
never observe isolated pentagons, heptagons or 5-7 dislocations. Contrary to graphite where 
melting is initiated by interplanar covalent bond formation, in graphene it seems that 5-7 
clusters act as nuclei for the melting. By close inspection, we find that regions with 5-7 
clusters favor the transformation of three hexagons into two pentagons and one octagon 
that we never see occurring in the regular hexagonal lattice far from the 5-7 clusters. The 
large bonding angle in octagons, in turn, leads to the proliferation of larger rings. Due to 
the weakening of the bonds with small angles in the pentagons around them, these larger 
rings tend to detach from the lattice and form chains. 

When melting is completed the carbon chains form an entangled 3D network with a sub- 
stantial amount of three-fold coordinated atoms, linking the chains. The molten phase is 
similar to the one found for fullerenes and nanotubes(see^ and references therein). There- 

16 



FIG. 8. Structure of graphene in the first phase of melting (top) and when molten (bottoni) at 
T = 5000 K. Adapted from RefP 

fore, the structure of the high temperature phase reminds rather a polymer gel than a simple 
liquid, a quite amazing fact for an elemental substance. 

The closest system to graphene is graphite. The melting temperature T m of graphite has 
been extensively studied experimentally at pressures around 10 GPa and the results present 
a large spread between 4000 K and 5000 K^. With LCBOPII, free energy calculations 
give T m = 4250 K, almost independent of pressure between 1 and 20 GPa (see references 
irP^l. At zero pressure, however, graphite sublimates before melting at 3000 Kp^. Monte 
Carlo simulations with LCBOPII at zero pressure show that, at 3000 K, graphite sublimates 
through detachment of the graphene layers. The melting of graphene in vacuum that we have 
studied here can be thought of as the last step in the thermal decomposition of graphite. 
Interestingly, formation of carbon chains has been observed in the melt zone of graphite 
under laser irradiatiorpS Although the temperature T = 4900 K of spontaneous melting 
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represents an upper limit for T m , our simulations suggest that T m of graphene at zero pressure 
is higher than that of graphite. 

IV. CONCLUSIONS 

We have shown by comparing the results of atomistic simulations to the theory of mem- 
branes based on a continuum approach, that graphene can indeed be considered a prototype 
of 2D membrane and that atomistic studies can be used to evaluate accurately the scaling 
properties, including scaling exponents and cross-over behavior. Conversely, the melting of 
graphene is determined rather by the peculiarities of the carbon-carbon bond and the high 
stability of carbon chains than as a generic model for melting in 2D. 
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